#ifndef __RGM_H__
#define __RGM_H__

#include <deque>
#include <vector>
#include <string>
#include <cassert>

#include "Image.H"
#include "Util.H"
#include "IO.H"

//putting namespaces in header files might not be the best solution
using namespace std;
using namespace pip;

// Arrays with offsets to 2D neighborhood.
// Taking the first 4 gives 4-neighborhood, taking all gives 8-neighborhood.
static const int dnbx2D[8] = {1, 0,-1, 0, 1, 1,-1,-1};
static const int dnby2D[8] = {0, 1, 0,-1, 1,-1, 1,-1};
// Arrays with offsets to 3D neighborhood.
// Taking the first 6 gives 6-neighborhood, taking first 18 gives 18-neighborhood,
// taking all gives 26 neighborhood.
static const int dnbx3D[26] = {1, 0, 0,-1, 0, 0,   1, 1, 1, 1, 0, 0,-1,-1,-1,-1, 0, 0,   1, 1, 1,-1,-1,-1, 1,-1};
static const int dnby3D[26] = {0, 1, 0, 0,-1, 0,   1,-1, 0, 0, 1, 1, 1,-1, 0, 0,-1,-1,   1,-1, 1, 1, 1,-1,-1,-1};
static const int dnbz3D[26] = {0, 0, 1, 0, 0,-1,   0, 0, 1,-1, 1,-1, 0, 0, 1,-1, 1,-1,   1, 1,-1, 1,-1, 1,-1,-1};

class Pixel
{
        private:
            int x, y, z;

        public:
        Pixel(int x = 0, int y = 0, int z = 0) {set(x,y,z);};

        template <class T>
        T get(const Image<T> &im) const {return im(x,y,z);}

        template <class T>
        void set(Image<T> &im,T val) const {im(x,y,z) = val;}

        void get(int &ix, int &iy, int &iz) const {ix=x;iy=y;iz=z;}

		bool negative()  {return (x < 0);}

        void set(int ix, int iy, int iz = 0) {x=ix;y=iy;z=iz;}

		Pixel& operator=(const Pixel &rhs) {if (this != &rhs) { x=rhs.x; y=rhs.y; z=rhs.z; } return *this; }

        template <class T>
        void getNeighbors(const Image<T> &im,vector<Pixel> &v, int n=8) const
        {
                v.clear();
                if(n==4||n==8){//2D
                        for(int i=0; i<n; ++i){
                                Pixel p = *this;
                                p.x+=dnbx2D[i];
                                if(p.x<0 || p.x>=im.dim(0)) continue;
                                p.y+=dnby2D[i];
                                if(p.y<0 || p.y>=im.dim(1)) continue;
                                //if we get here, the position is in the image
                                v.push_back(p);
                        }
                }
                else if(n==6||n==18||n==26){//3D
                        for(int i=0; i<n;++i){
                                Pixel p =*this;
                                p.x+=dnbx3D[i];
                                if(p.x<0 || p.x>=im.dim(0)) continue;
                                p.y+=dnby3D[i];
                                if(p.y<0 || p.y>=im.dim(1)) continue;
                                p.z+=dnbz3D[i];
                                if(p.z<0 || p.z>=im.dim(2)) continue;
                                v.push_back(p);
                        }
                }
                else
                assert(false);//illegal value of n
        }

};

template <class T, class U>
int connect(
const Image<T> &in,
Image<U> &out,
U target, // value for connected object in out
int n, // connectivity
int x, int y,
int z = 0 // optional
)
{
        //if out is not the same size as in,
        //resize it and fill with 0
        if (out.dimension()!=in.dimension()) {
                out.resize(in.dimension());
                out.fill(U(0));
        }

        Pixel seed(x,y,z);
        T val = seed.get(in);
        seed.set(out,target);
        deque<Pixel> q;
        q.push_back(seed);
        int count = 1;//first pixel has been added to q already
        while (!q.empty()) {
                deque<Pixel>::iterator i = q.begin();
                vector<Pixel> nn;
                i->getNeighbors(in,n,nn);
                for (vector<Pixel>::iterator j=nn.begin();j!=nn.end();++j) {
                        if ((j->get(in)==val)&&(j->get(out)!=target)) {
                                q.push_back(*j);
                                j->set(out,target);
                                ++count;
                        }
                }
                q.pop_front();
        }

        return count;
}

template <class T>
void connectedComponentLabeling(
const Image<T> &in, // input image
Image<float> &out, // output image
int n, // connectivity
vector<int> &sizes) // returns the sizes of the objects
{
        //check sizes of output image
        if (out.dimension()!=in.dimension()) {
                out.resize(in.dimension());
                out.fill(0);
        }

        sizes.clear();
        float label = 1;//label for connected objects

        for(int iz=0; iz<in.dim(2); ++iz){
                for(int iy=0; iy<in.dim(1); ++iy){
                        for(int ix=0; ix<in.dim(0); ++ix){
                                if ((in(ix,iy,iz)!=T(0))&&(out(ix,iy,iz)==0)){
                                        //pixel is not background and not in already grown object
                                        int count = connect(in,out,label,n,ix,iy,iz);
                                        sizes.push_back(count);
                                        ++label;
                                }
                        }
                }
        }
}

template <class T, class U, class GrowCheck>
int grow(
const Image<T> &in,
Image<U> &out,
GrowCheck gc,
U target, // value for connected object in out
int n, // connectivity
int x, int y,
int z = 0)
{
        //resize out, if we must
        if (out.dimension()!=in.dimension()) {
                out.resize(in.dimension());
                out.fill(U(0));
        }
        Pixel seed(x,y,z);
        seed.set(out,target);
        vector<Pixel> nn;
        deque<Pixel> q;
        q.push_back(seed);
        int count = 1; //first pixel has been added
        while (!q.empty()) {
                deque<Pixel>::iterator i = q.begin();
                i->getNeighbors(in,n,nn);
                for (vector<Pixel>::iterator j=nn.begin();j!=nn.end();++j) {
                        if (gc(j->get(in))&&(j->get(out)!=target)) {
                                q.push_back(*j);
                                j->set(out,target);
                                ++count;
                        }
                }
                q.pop_front();
        }
        return count;
};

template <class T>
class ConnectCheck
{
        private:
        T t;
        public:
        ConnectCheck(const T &val) : t(val) {}
        bool operator()(const T &val) {return val==t;}
};

template <class T, class U>
int connect2(
const Image<T> &in,
Image<U> &out,
U target, // value for connected object in out
int n, // connectivity
int x, int y,
int z = 0)
{
        ConnectCheck<T> gc(in(x,y,z));
        return grow(in,out,gc,target,n,x,y,z);
}

template <class T>
class RegionGrowUpCheck
{
        private:
        T t;
        public:
        RegionGrowUpCheck(const T &val) : t(val) {}
        bool operator()(const T &val) {return val<=t;}
};

template <class T>
class RegionGrowDownCheck
{
        private:
        T t;
        public:
        RegionGrowDownCheck(const T &val) : t(val) {}
        bool operator()(const T &val) {return val>=t;}
};

template <class T, class U>
int regionGrowUp(
const Image<T> &in,
Image<U> &out,
T threshold,
U target,
int n, // connectivity
int x, int y,
int z = 0)
{
        RegionGrowUpCheck<T> rgu(threshold);
        return grow(in,out,rgu,target,n,x,y,z);
}

template <class T, class U>
int regionGrowDown(
const Image<T> &in,
Image<U> &out,
T threshold,
U target,
int n, // connectivity
int x, int y,
int z = 0)
{
        RegionGrowDownCheck<T> rgd(threshold);
        return grow(in,out,rgd,target,n,x,y,z);
}

template <class T, class U>
int regionGrow(const Image<T> &in,
Image<U> &out,
T threshold,
U target, // value for connected object in out
int n, // connectivity
int x, int y,
int z = 0)
{
        T val = in(x,y,z);
        if (val<=threshold) {
                cout << "Value at start point is " << val << ", Threshold is " << threshold
                << " so everything below or equal to the threshold will be grown." << endl;
                return regionGrowUp(in,out,threshold,target,n,x,y,z);
        }
        else {
                cout << "Value at start point is " << val << ", Threshold is " << threshold
                << " so everything above or equal to the threshold will be grown." << endl;
                return regionGrowDown(in,out,threshold,target,n,x,y,z);
        }
}

template <class T> void hysteresisThresholding(const Image<T> &in, // input image
Image<unsigned char> &out, // output image
int n, // connectivity
T tHigh,
T tLow)
{
        //check sizes of output image
        if (out.dimension()!=in.dimension())
        out.resize(in.dimension());
        out.fill(0);

        for(int iz=0; iz<in.dim(2); ++iz)
        for(int iy=0; iy<in.dim(1); ++iy)
        for(int ix=0; ix<in.dim(0); ++ix)
        if ((in(ix,iy,iz)>=tHigh)&&(out(ix,iy,iz)==0))
        regionGrowDown(in,out,tLow,(unsigned char) 1,n,ix,iy,iz);
}

#endif